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ABSTRACT 

We use a high resolution ACDM numerical simulation to calculate the mass function 
of dark matter haloes down to the scale of dwarf galaxies, back to a redshift of fifteen, 
in a 50 ft. _1 Mpc volume containing 80 million particles. Our low redshift results allow 
us to probe low a density fluctuations significantly beyond the range of previous 
cosmological simulations. The Sheth and Tormen mass function provides an excellent 
match to all of our data except for redshifts of ten and higher, where it overpredicts 
halo numbers increasingly with redshift, reaching roughly 50 percent for the 10 10 — 
1O U M0 haloes sampled at redshift 15. Our results confirm previous findings that the 
simulated halo mass function can be described solely by the variance of the mass 
distribution, and thus has no explicit redshift dependence. We provide an empirical fit 
to our data that corrects for the overprediction of extremely rare objects by the Sheth 
and Tormen mass function. This overprediction has implications for studies that use 
the number densities of similarly rare objects as cosmological probes. For example, 
the number density of high redshift (z ~ 6) QSOs, which are thought to be hosted by 
haloes at 5a peaks in the fluctuation field, are likely to be overpredicted by at least a 
factor of 50%. We test the sensitivity of our results to force accuracy, starting redshift, 
and halo finding algorithm. 

Key words: galaxies: haloes - galaxies: formation - galaxies: clustering - cosmology: 
theory - cosmology: dark matter 



1 INTRODUCTION 

Cold dark matter models with a cosmological constant 
(ACDM) are successful in explaining a wide array of kine- 
matic and structural properties of the observed universe. A 
critical test of the ACDM model is how well it predicts the 
abundance of dark matter haloes, which serve as hosts for 
observable clusters, groups, and galaxies. Simulations that 
resolve haloes out to high redshift can be used to model 
the evolution of the numbers of observable high redshift 
objects, their progenitors, and their evolved descendants. 
Lyman break galaxies, for example, observed at redshifts 
out to z ~ 4 (e.g. Steidel et al. 1996) are likely progeni- 
tors of groups or clusters (Governato et al. 1998; Governato 
et al. 2001), and we are able to model their numbers over 
their entire observable lifespan. Many objects, however, lie 
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outside the realm that can presently be simulated because 
their number densities, masses, or redshifts are too extreme. 
For example, to model a reasonably sized sample of the 
hosts of the highest redshift (z ~ 6) quasi-stellar objects 
(QSOs), whose number density has been measured by the 
Sloan Digital Sky Survey (Fan et al. 2001), would require 
a simulation with volume roughly as large as the observ- 
able universe with (so far) prohibitively high particle num- 
bers. Simulations of adequate resolution and volume could 
be used, in principle, to estimate the host masses of such 
rare QSOs, or to estimate cosmological parameters after as- 
suming host masses, by matching predicted and observed 
number densities. However, with much less computational 
effort, by modeling smaller cosmological volumes at higher 
redshift, we are able to test analytic mass functions in the 
same regime of rare density enhancements. Generally, rare 
density peaks correspond to high values of M/M*, where 
M* is the 'characteristic' mass of a typical collapsing halo 
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at that epoch (to be discussed later). By modeling haloes 
over a wide range of M/M*, we can constrain analytic mass 
functions for ranges of redshift and mass that have not yet 
been simulated. This is possible because analytic mass func- 
tions are generally derived from assumptions of how linear 
density fluctuations lead to halo collapse, and thereby have 
no explicit redshift dependence. High redshift QSO hosts, 
galaxy progenitors, and perhaps even the first generation of 
stars are all examples of high M/M* objects whose mass 
functions can presently only be calculated analytically. 

Press & Schechter (P-S, 1974) developed an analytical 
framework that predicts the number and formation epoch of 
dark matter haloes. In P-S theory, as the universe evolves, 
linear density fluctuations grow gravitationally until they 
reach a critical spherical overdensity, at which time non- 
linear gravitational collapse is assumed to occur. Cosmolog- 
ical numerical simulations have shown the P-S framework 
to be approximately correct, but P-S theory consistently 
underpredicts the number of high mass haloes and overes- 
timates the number of haloes less than about M, (e.g. Efs- 
tathiou et al. 1988; Gross et al. 1998; Governato et al. 1999; 
Jenkins et al. 2001; White 2002) even when merging of dark 
matter haloes is included in predictions (Bond et al. 1991; 
Bower 1991; Lacey & Cole 1993; Gardner 2001), though the 
high mass end fits well if the finite size of haloes is taken 
into account (Yano, Nagashima, & Gouda 1996; Nagashima 
2001). Ellipsoidal halo collapse models (e.g. Monaco 1997ab; 
Lee & Shandarin 1998; Sheth, Mo, & Tormen 2001) yield 
much more robust predictions than the conventional spheri- 
cal collapse models, and are in excellent agreement with em- 
pirical fits by Sheth & Tormen (1999). Monaco et al. (2001), 
using the semi-analytic code PINOCCHIO, which uses a per- 
turbative approach, show that the dark matter halo distri- 
bution can be accurately predicted at much lower computa- 
tional cost, on a point-by-point basis, from a numerical real- 
ization of an initial density field. Jenkins et al. (2001) utilize 
a large set of simulations of a range of volumes and cosmolo- 
gies (Jenkins et al. 1998; Governato et al. 1999; Evrard et 
al. 2002) to test the Sheth & Tormen (S-T) mass function 
over more than four orders of magnitude in mass, and out to 
a redshift of 5, finding good agreement with the S-T function 
down to their resolution limit of ~ 3 x 10 11 Mq , except for an 
overprediction by the S-T function for haloes at rare density 
enhancements. In this study, we probe previously untested 
regimes of the mass function by simulating a volume that 
resolves haloes down to the scale of 10 10 Mq dwarfs, in a 
cosmological environment, allowing us to sample the mass 
function back to z~15. Our paper is outlined as follows: in 
section 2, we describe our simulations and numerical tech- 
niques; in section 3, we review the analytic theory; then we 
discuss our results and compare with previous work in sec- 
tion 4; we conclude with a discussion of the implications of 
our work. 



2 THE SIMULATIONS 

Our cosmology is the presently favored ACDM with A = 0.7 
and Q m = 0.3. We use the parallel tree gravity solver PKD- 
GRAV (Stadel 2001) to simulate 81 x 10 6 (432 3 ) dark matter 
particles from a starting redshift, zo, of 69. We then resim- 
ulate the same volume but with zo =139, and evolve this 



volume to z=7; which allows us to consider results at higher 
redshift than with the zo =69 run. In order to simulate the 
highest possible mass resolution, we employ a volume of 50 
/i -1 Mpc on a side. Our particle mass is 1.3 x lO 8 ft _1 M0 
allowing us to resolve haloes down to less than 10 10 /i _1 Mq 
with 75 particles. Our force resolution is 5 /i _1 kpc. We use a 
cell opening angle of O < 0.8 at low redshift, and O < 0.7 at 
z > 2. In section 4.2, we discuss tests that confirm that our 
choices of initial redshift, softening, and high redshift open- 
ing angle are adequate (see Table 1). We use a "multistep- 
ping" approach, where particles in the highest density re- 
gions underg o 16 ,000 timesteps. Timesteps were constrained 
to 5t < 0.2a/ e/a, where e is the softening length and a is the 
magnitude of the acceleration of a given particle. We normal- 
ize the density power spectrum of our initial conditions such 
that as, the rms density fluctuation of spheres of 8 /i _1 Mpc 
extrapolated to redshift of zero is 1.0, consistent with both 
the cluster abundance (see e.g. Eke, Cole, & Frenk 1996 and 
references therein) and the COBE normalization (e.g. Ratra 
et al. 1997). To set our initial conditions, we use the Bardeen 
et al. (1986) transfer function with 7 = Q m ah, where h is 
the hubble constant in units of 100 kms~ Mpc—1. 

2.1 Halo Identification 

In order to identify haloes in our simulation, we use both 
the friends-of-friends (FOF) algorithm (Davis et al. 1985), 
and the spherical overdensity (SO) algorithm (Lacey & Cole 
1994). The FOF halo finder uses a "linking length" , 11, to link 
together all neighboring particles with spacing closer than 
II as members of a halo. SO identifies haloes by identifying 
spherical regions with the expected spherical overdensities 
of virialized haloes. For our SO haloes, we first use SKID 
(Stadel 2001) to identify all bound haloes, including those 
that are subhaloes within larger haloes. Next, we grow a 
sphere outward from each SKID center until it just con- 
tains the virialized overdensity. Finally, we iterate by grow- 
ing spheres outward from all neighboring SKID haloes until 
we have identified the center of mass for each SO halo. For 
our SO criterion, we use the virial overdensity predicted by 
the spherical collapse tophat model of Eke, Cole, & Frenk 
(1996). In a ACDM universe, the virial overdensity, in units 
of critical density, declines from an asymptotic value of 178 
as cosmic time increases and f2 m drops below 1; for red- 
shift of zero, the ACDM overdensity A V i r is 100 (Kitayama 
& Suto 1996). We exclude haloes that contain less than 64 
particles, which is conservative, as it is more than the esti- 
mated 20 or 30 particles needed for a robust halo identifica- 
tion based on resolution tests (Jenkins et al. 2001; Governato 
et al. 1999). 

We utilize the FOF algorithm for the bulk of our anal- 
yses since it is reasonably robust and computationally effi- 
cient. Our FOF // choice is 0.2 for all redshifts (except when 
matching 11 from previous studies); this approach has been 
shown to be sound for a range of cosmologies by Jenkins 
et al. (2001), although the evolution of A V i r in the spheri- 
cal collapse tophat model implies that 11 should range from 
=0.164 at z=0 to II =0.2 at high redshift in ACDM cos- 
mology (Lacey & Cole 1994; Eke, Cole, & Frenk 1996; Jenk- 
ins et al. 2001). To test the sensitivity of the mass function 
to halo selection criteria, we apply both FOF 11 =0.164 and 
SO to our data at various redshifts, discussed in section 4.1. 



© 2003 RAS, MNRAS 000, [Dig] 



Evolution of the Mass Function 3 



3 ANALYTIC THEORY 

Analytic P-S formalism yields the following (See Jenkins et 
al. 2001, whose notation we adopt in this section, and ref- 
erences therein for further discussion, which we summarize 
here) : 



/(a; P-S) = y exp 

" 7r a 



2 8 C 



2a 2 



(1) 



where S c is the threshold spherical linear overdensity above 
which a region will collapse. 8 C depends only weakly on the 
cosmological parameters and redshift (e.g. More, Heavens, 
& Peacock 1986; Jenkins et al. 2001), and S c = 1.686 for 
Q.0 = 1. In our analysis, we assume that 5 C = 1.686 at all 
redshifts. a 2 (A4, z) is the variance of the linear density field, 
smoothed with a spherical top-hat filter enclosing mass M, 
and is calculated from the linear density power spectrum 
P(k), extrapolated to z = 0: 



k 2 P(k)W 2 (k;M)dk, 



(2) 



where W(k; M) is the Fourier-space top-hat filter, and b(z) 
is the growth factor of linear perturbations normalized to 
unity at z=0 (Peeble 1993). In the P-S formalism, all mass 
is contained in haloes: 



/(a;P-S)dlncr" 



(3) 



The mass function f(a, z) can be related to the number den- 
sity, n(M, z), of haloes with mass less than M: 



M dn(M,z) 
po(z) diner -1 



(4) 



where po(z) is the mean density of the universe at that time. 

The S-T model is a modification to the P-S model based 
on empirical fits to simulations (Sheth & Tormen 1999), has 
been shown to reproduce simulation results substantially 
better than P-S (e.g. Jenkins et al. 2001; White 2002), and 
is theoretically justified in that it matches P-S formalism 
derived with ellipsoidal halo collapse models (Sheth, Mo, & 
Tormen 2001): 



/(a; S-T) = A 



r 



• exp 



aS 2 
2a 2 



(5) 



where A=0.3222, a = 0.707 and p = 0.3. Jenkins et 
al. (2001) offer an empirical fit using high resolution sim- 
ulations of a range of cosmologies. Their fit is constructed 
in the / — ln(<r -1 ) plane, which has the advantage of being 
invariant with redshift: 

/(In a' 1 ) = 0.315 exp [ - | In a' 1 + 0.61| 3 ' 8 ] . (6) 

The Jenkins et al. function adjusts for an overprediction by 
the S-T function for the rare objects at large lner -1 , and is 
calibrated for the range — 1.2 < In a -1 < 1.05, which corre- 
sponds to masses down to approximately 3x10 h~ Mq at 
present epoch, and includes haloes out to z=5, with FOF 
fixed at 11 =0.2. 

In each of these analytic functions, virialized haloes 
have a characteristic mass, M*(z): 




a{M*{z)) = S c 



(7) 



Figure 1. The curves correspond to the halo mass of ncr fluctu- 
ations in the density field, given by cr(M, z) = S c /n, with n =1, 
2, 3, 4, and 5, from bottom to top. The characteristic mass, M*, 
given by a(M„(z)) = <5 C , is the la curve. The area between the 
long dashed lines is sampled with poisson errors of less than 20% 
in our data. 



and a(M, z) — o(M, z = 0)b(z). b(z) evolves as (1 + z)~ in 
an fio = 1 universe, and more slowly in a ACDM universe. 
a(M) decreases slowly with increasing mass, which leads to 
the steep redshift dependence of M*(z), shown in Fig. 1 for 
the ACDM universe, and results in a broad mass spectrum of 
collapsed objects. An important test of analytic mass func- 
tions is the accuracy of their predictions at high values of 
M/M* and z. At low redshift, high M/M, haloes would have 
unrealistically large masses, but at high redshift the steep 
evolution of M* puts high M/M, objects well into the realm 
of simulations. Furthermore, the evolution of halo masses 
that correspond to high a density enhancements means that 
rare haloes are most easily simulated at high redshift (Fig. 
1). Our simulations model the mass function of haloes lying 
at up to 4er density fluctuations. 



4 EVOLUTION OF THE MASS FUNCTION 

In Fig. 2, we compare the analytic version of the P-S, the 
S-T, and the Jenkins et al. mass function with our sim- 
ulation results at several redshifts. The S-T function pro- 
vides the best fit to our simulation with excellent agreement 
at all masses and redshifts except for our highest redshift 
outputs. The P-S function overpredicts substantially every- 
where except at the low and high mass extremes. And the 
Jenkins et al. mass function fits much of our data well at 
z=0, but diverges from our simulation results once well be- 
low the limit of its empirical fit of ln<j _1 = —1.2, which 
corresponds to ~ 4 x 10 11 /i" 1 Mq with as = 1.0. Note that 
a(M,z) is sensitive to as, so for a as = 0.9 ACDM model, 
which was used in many of the simulations that were part 
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Figure 2. Comparison of the mass function per decade of mass. 
Data points are our ACDM simulation results with lcr poisson 
error bars. Throughout the paper, haloes are identified using a 
FOF 11 =0.2, unless otherwise specified; only haloes with at least 
64 particles are considered. In our plots, we plot the median halo 
mass in a bin. Solid curves are the Sheth & Tormen function at 
z=0, 2, 5, 8, & 15. Short dashed curve is the Jenkins et al. "uni- 
versal" mass function (Eqn. 6), which diverges when extrapolated 
well below its original lower mass limit of ~ 4 X l0 11 h,- 1 M e for 
our as = 1.0 model (see text). Long dashed curve is the Press &; 
Schechter function. 



of the Jenkins et al. fit, lna 1 = —1.2 would correspond to 
only ~ 2 x 10 11 r 1 M e . 

In Fig. 3, we show the evolution of the mass function 
over all of our redshifts compared to its S-T predicted evolu- 
tion. The S-T function provides an excellent fit to our data, 
except at very high redshifts, where it significantly over- 
predicts the halo abundance. At all redshifts up to z=10, 
the difference is ;$ 10% for each of our well sampled mass 
bins. However, the S-T function begins to overpredict the 
number of haloes increasingly with redshift for z £ 10, up 
to ~ 50% by z=15. The simulation mass functions appear 
to be generally steeper than the S-T function, especially at 
high redshifts. In Fig. 4, we show the evolution of the mass 
function over all of our redshifts as a function of M/K'h. This 
highlights the remarkable accuracy of the S-T mass function 
over more than 10 decades of M/M*. 

In Fig. 5, we plot the mass function for all of our out- 
puts in the / — ln(cr ) plane. Large values of ln<r~ cor- 
respond to rare haloes of high redshift and/or high mass, 
while small values of lntr^ 1 describe haloes of low mass and 
redshift combinations. In Fig. 6, we compare our simulated 
mass function with the S-T prediction, by plotting the resid- 
uals over our entire range. We limit plotted data to bins with 
poisson errors of less than 20%. Remarkably, the S-T func- 
tion fits our simulated mass function to better than 10% 
over the range of -1.7 ^ lna -1 ^ 0.5. We are unaware of 
any previous studies that probe the mass function down to 
such low values of In <r _1 in a cosmological environment. The 
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Figure 3. Fractional difference between our simulated mass func- 
tion (FOF 11 =0.2) and the S-T prediction. Poisson error bars 
shown. 
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Figure 4. Number density vs. M/M*. Data points with poisson 
error bars are our simulation results. Curves are S-T predictions. 
Redshifts plotted are (from left to right) 0, 1., 2., 3., 4., 5., 6.2, 
7.8, 10., 12.1, 14.5 



S-T function appears to significantly overpredict haloes for 
lna -1 ^ 0.5. This is the same overprediction seen in the 
number density for z ^ 10 in figs. 3 and 4. The large ap- 
parent scatter of the mass function for lna- -1 ^ 0.5 is due 
to larger poisson errors in this range. For large lncr -1 , we 
estimate the uncertainty in the mass function due to cosmic 
variance by estimating the contribution of linear fluctua- 
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tions on the scale of the box size. Cosmic variance can have 
large effects on results since it has the potential to increase 
or decrease the mass function in multiple mass bins simulta- 
neously, and cosmic variance is difficult to quantify without 
a large set of simulated volumes. In our estimation, we use 
a second order Taylor expansion of f(a; S— T), ignoring the 
first order term since we have imposed the average density 
of our simulation to be Q m . The resulting estimate for the 
uncertainty in the mass function due to cosmic variance is 

A/, c .„. ~ 9 Q0>2~ ^cltboxi which we evaluate numerically. 
In our well sampled, low redshift mass bins, A/ lC .„. is negli- 
gible. For our highest redshift results, A/ iC .„. is smaller than 
our poisson error limit of 20% for figures 6 and 7, even for 
our highest mass bins. However, Af tC , v , approaches our pois- 
son error limit of 20% for our z =14.5 output. For our z =10 
and z =12 outputs, A/ j0 .„. is less than 10% in bins where 
poisson errors are less than 10%, which is the case for most 
of our bins at that redshift. Thus, while cosmic variance is a 
significant source of error where the mass function is steep- 
est, it is unlikely to entirely account for our discrepancy with 
the S-T function. We note that several of our z < 2 points 
lie roughly 3<r above the mean; this is actually just one mass 
bin plotted repeatedly at different redshifts, and so is not en- 
tirely surprising. By careful examination of ranges of lncr -1 
where outputs of different redshifts overlap, we verify that 
the magnitude of the S-T overprediction at high values of 
lncr" 1 is consistent with being a function purely of lncr - 
rather than redshift, a natural consequence of the fact that 
the mass function is self similar in time (e.g. Efstathiou et 
al. 1988; Lacey & Cole 1994; Jenkins et al. 2001). Jenkins 
et al. (2001) also find an overprediction by the S-T func- 
tion for lncr" 1 £ 0.75, which with their larger simulation 
volumes, corresponded primarily to objects of z^2 and of 
much higher mass. Additionally, Jenkins et al. find the mass 
function to be invariant with redshift within their own re- 
sults. In Fig. 7, we compare a subset of our data with the 
Jenkins et al. "universal" mass function. We note that when 
extrapolated to lncr" 1 ^ —1-4, well below its empirical fit 
range of —1.2 ^ lncr" 1 ^ 1.05, the Jenkins et al. func- 
tion diverges from our results, reflecting the fact that it is 
of a form not ideally suited for extrapolation. Where our 
data overlap (lncr" 1 ^ —1.2), we find generally good agree- 
ment, although we have ~ 20% fewer haloes and a some- 
what steeper mass function at lncr" 1 0.25. Over the range 
of -1.4 lncr" 1 0.75, the Jenkins et al. mass function 
matches our data to within ~ 20%. 

We consider the possibility that this difference between 
our data and the Jenkins et al. fit could be due to differences 
in the effective slope of the power spectrum, n e //, where 
P(k) oc k n "ff . Applying a simple power law to Eqn. 2 yields 
cr 2 cx M~( Ue ff +3 )/ 3 , which can be reparameterized as 



6 



diner 
'din A/ 



(8) 



(Jenkins et al. 2001). Fig. 8 shows the n e ff vs. lncr -1 pa- 
rameter space that our results cover. Our data generally has 
a steeper n e ff than Jenkins et al., though there is some 
overlap, especially at lower redshifts. Since the slope of the 
linear power spectrum is invariant with redshift for a given 
k, n e ff(z) is constant for a given mass, meaning that for our 
lowest mass haloes we sample nearly all of our /(cr) at an 
n e ff ~ —2.3. If we consider n e //(cr = 0.5), where the power 




Figure 5. Mass function plotted in redshift independent form for 
all of our outputs as in Fig. 4. Solid curve is S-T prediction. 



spectrum begins to go nonlinear, then we find that our re- 
sults differ significantly from the Jenkins et al. function only 
where n e ff at the nonlinear scale is steepest. In particular, 
the disagreement is worst at z^lO where n e //(cr = 0.5) ex- 
ceeds their maximum value of -2.26. Thus, there is a pos- 
sibility that our steeper values of n e // results in an f(a) 
with a slightly different shape, which might account for the 
difference with prior results, but a larger set of simulations 
with still steeper n e ff would be needed to clearly show any 
such dependence. 

Since no mass function that we have considered is ac- 
curate for the entire range of our data, we consider the pos- 
sibility of an empirical adjustment to the S-T function. We 
insert a crude multiplicative factor to the S-T function as 
follows, with S c = 1.686 and FOF 11 =0.2 (Fig. 6): 



/(cr) = /(a; S-T) 



ea:p[-0.7/(cr[cosh(2cr)] 5 )] 



(9) 



valid over the range of -1.7 ^ lncr" 1 ^ 0.9. The result- 
ing function is virtually identical to the S-T function for 
all — oo ^ lncr 0.4. At higher values of lncr -1 , this func- 
tion declines relative to the S-T function, reflecting an un- 
derabundance of haloes that becomes greater with increas- 
ing lncr -1 . For -1.7 ^ lncr" 1 ^ 0.5, Eqn. 9 matches our 
data to better than 10% for well-sampled bins, while for 0.5 
^ lncr -1 ^ 0.9, where poisson errors are larger, our data is 
matched to roughly 20%. We must caution that though Eqn. 
9 is a good fit to our data, it differs from the S-T function 
in the regime where poisson and cosmic errors are highest, 
and where our results are most prone to potential numerical 
errors because of the steepness of the mass function. Our 
results are more robust in the low lncr -1 regime. Note that 
in the Eqn. 9 fit, not all mass belongs to a halo, so Eqn. 3 
is not valid. 
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Figure 6. Residuals between S-T prediction and our results for 
the mass function of Fig. 5 (with FOF 11=0.2). Solid straight line 
is the S-T function. Dashed line is our empirical adjustment to 
the S-T function. 



Figure 8. Parameter range covered by our results. Each mass 
bin for which we have poisson statistics of better than 20% error 
is shown. The vertical line denotes a = 0.5. 



< 




Figure 7. Residuals between Jenkins et al. mass function and 
our results for the mass function of Fig. 5. Arrows encompass the 
range of data used in the Jenkins et al. empirical fit, which is 
denoted by the solid straight line. 



4.1 Friends-of-Friends (FOF) versus Spherical 
Overdensity (SO) 

Other authors have noted the advantages and disadvantages 
of the FOF and SO algorithms (see Jenkins et al. 2001 and 
references therein) . The FOF method has the advantage that 
it can identify haloes of any shape as long as their minimum 



local number density is at least roughly 1/b 3 , and FOF is 
generally computationally cheaper than SO. However, FOF 
can sometimes spuriously link together haloes that lie close 
together within a filament (see e.g. Governato et al. 1997). 
In Fig. 9, we compare FOF mass functions of our simulation 
with the corresponding SO mass functions for a range of red- 
shifts. Note that the vertical axes are somewhat arbitrary 
for the z=10 and z=15 outputs as these were made from the 
zo = 69 version of the simulation which had a somewhat sup- 
pressed mass function, which we discuss in section 4.2. The 
halo finders have excellent agreement at low redshifts, with 
differences of 10% over the range where the mass func- 
tion is well sampled. Differences in the high mass bins are 
due to a combination of different mass calculations for indi- 
vidual selected clusters as well as offset mass bins. The steep 
dropoff in the SO mass function for low masses is due to a 
our exclusion of SKID haloes of less than 64 particles as po- 
tential SO centers. At high redshifts, FOF 11 =0.2 produces 
a substantially higher mass function than FOF =0.164 
or SO, which are similar to each other, implying sensitivity 
to 11 for large In cr _1 , probably because the mass function 
is steep there, and thus sensitive to halo selection criteria. 
To verify that the discrepancy of FOF U =0.2 with SO is 
not due simply to our choice of using skid haloes as our ini- 
tial SO centers, we have included a high redshift (z=10) SO 
mass function which uses FOF haloes as initial SO centers; 
our choice of centers from which to grow our SO spheres has 
no effect on the SO mass function. A visual inspection in 
which halo members are "marked" , reveals that at high red- 
shift, FOF 11 =0.2 links together some neighboring haloes 
connected by filaments, but FOF also identifies some in- 
dividual haloes (often of highly elongated shape) that are 
missed by SO, so neither algorithm is ideal. The overpredic- 
tion of the S-T function for rare objects worsens somewhat 
if we use SO derived mass functions. Using a linking length 
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Figure 9. Comparison of SO and FOF mass functions from our 
simulation. Filled squares (connected by dashed lines) are for FOF 
11 =0.2 haloes; filled triangles (connected by dotted lines) are FOF 
11 =0.164 haloes. Open squares (solid lines) are SO haloes made 
with SKID centers. X's at z=10 are SO haloes made using FOF 
centers. SO overdensity criterion is from Kitayama & Suto 1996. 
Note that the z = 10 and z = 15 data lie too far below the S-T 
function in this plot because it is from our zo =69 run rather than 
our zo =139 run, which we use for the rest of the high redshift 
results of this paper. 

that varies with redshift in an attempt to match the varying 
overdensity of virialized haloes, would have little effect on 
our mass function, since at low redshift the mass function 
is insensitive to U, and at high redshift tt m — 1, implying 
11 =0.2 (Davis et al. 1985; Lacel & Cole 1994). However, 
had we sampled large In cr _1 at low redshift, adjusting 11 
to match virial overdensity would likely have a significant 
effect. 



4.2 Numerical Tests 

We check that the disagreement of S-T which appears at 
high redshift is not a result of delayed halo collapse due to 
numerical errors that can be caused by too large of a gravita- 
tional softening length. We make additional checks address- 
ing potential numerical errors caused by mapping particles 
with Zel'dovich displacements (Zel'dovich 1970) onto a par- 
ticle grid; an insufficiently high starting redshift could delay 
collapse of the first haloes (e.g. Jenkins et al. 2001). If ini- 
tial conditions are set with some regions having overdensities 
high enough to already be in the non-linear regime, then the 
linear Zel'dovich mapping can not account for shell- crossing 
wherein mass piles up as it flows toward overdensities. The 
effects of either of these error sources, if present, should have 
evolved away by lower redshifts, since the tiny fraction of 
matter that is in haloes at such high redshifts is soon incor- 
porated into clusters or large groups. 

Table 1 lists our test runs, each of which consists of an 
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Figure 10. The high redshift mass functions for runs with soften- 
ing lengths of 5 h _1 kpc and 2.5 h~ x kpc. Both runs were started 
from a redshift of 69. 

identical 432 3 particle volume with identical random waves, 
and is evolved to z=7. Fig. 10 shows the mass function for 
z~7-15 for our low softening test run, started from zo =69, 
and plotted relative to the S-T function along with our 
zo =69, 5 /i -1 kpc data. Halving the softening to 2.5 ft _1 kpc 
has no effect on the mass function. Fig. 11 shows the mass 
function for z~7-15 for our initial redshift test runs, all with 
5 /i _1 kpc softening. Lowering the initial redshift to 39 sub- 
stantially reduces the number of high redshift haloes, so 
we did not evolve the zo =39 run to z<10. The zo =139 
run matches the zo =279 run, indicating convergence, but 
the zo =69 run has a reduced mass function relative to 
the zo =139 run at redshifts z £ 12. By z=7, however, the 
zo =139 and zo =69 mass functions have converged, showing 
that evolving the simulation over an expansion factor of ~10 
from initial conditions is sufficient for mass function mea- 
surements. We consequently derive our low redshift (z<7) 
results throughout this paper from the zo =69, 5 /i _1 kpc 
simulation, and utilize the zo =139 run for our z^7 results. 
We make an additional test of zjs7 cell opening angle, which 
is used to determine how accurately long range gravitational 
forces are to be approximated. If the cell opening angle is 
too large, then artificial net forces will be incurred upon 
particles, which could cause "spurious" haloes to form. This 
effect is most likely to occur at high redshifts when grav- 
itational perturbations are small and force errors are frac- 
tionally larger. Decreasing the opening angle from 0=0.7 to 
O=0.5, for our zo =139 case, had no appreciable effect on the 
mass function for z^7. We test that the number of replicas 
used for our periodic boundaries, n r =1, is adequate. With 
zo =139, increasing n r from 1 to 2 has no effect on the mass 
function. Additionally, to test that our box size is adequate, 
we have verified that our z=0 mass function agrees with the 
mass function from larger, lower resolution volumes where 
they overlap (not included in Table 1). 
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Table 1. N-body simulation parameters, including test runs 

ZO r sa ft(A _1 kpc) 6(z>7) rVepiica Zei-oi^ed: 6(2<z<7) 0(z<2) 

Applied to z<7 results; SO vs. FOF test 
Applied to z^7 results 
Test 
Test 
Test 
Test 



69 

139 

139 

69 

39 

279 



5.0 
5.0 
5.0 
2.5 
5.0 
5.0 



0.7 
0.5 
0.5 
0.7 
0.7 
0.5 
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10 
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o.s 




Log lo (M/h-'M ) 

Figure 11. The high redshift mass functions for initial redshifts 
of 39, 69, 139, and 279. Each run has a softening of 5 /i _1 kpc. 
The zq =39 run was stopped at z =10. 



5 CONCLUSIONS 

Our results extend to lower masses and higher redshifts than 
the original empirical fit of the S-T mass function. The range 
of masses and redshifts over which the S-T mass function re- 
mains valid is quite remarkable, and though it does begin 
to break down at redshift £ 10 in our results, no other 
function matches its range and accuracy. It is not well un- 
derstood why the mass function can be described so well by 
solely a(M). Lacey & Cole 1994, in simulations with scale 
free power spectra, found evidence that the mass function 
depends on n e ff as we U as °", though they tested a wider 
range in n e // than in more recent CDM simulations. The 
generally good agreement of our data with previous work 
at much shallower n e // implies that any such dependence 
is weak, though it may be manifested in our results where 
we differ from the S-T function. As n e ff approaches -3, the 
growth of M, with time diverges, so any dependence of f(<r) 
on n e / / is most likely to occur near that regime. Simulations 
with finer mass resolution will be able to test yet steeper val- 
ues of n e ff. Based on the apparent trend of the S-T function 
to overpredict halo numbers for objects of greater rarity, we 
expect that the overprediction of the S-T function may con- 
tinue to worsen when even higher values of In a" 1 , or equiv- 



alently, when more extreme values of M/M, and redshift, 
are analyzed. Theoretical work that focuses on halo collapse 
in this high In a~ x range is needed to produce more robust 
predictions. 

Because the form of the mass function for low mass, 
low redshift haloes closely resembles a power law, we are 
cautiously optimistic that the S-T mass function will con- 
tinue to provide a good match to simulation data as lower 
values of In a -1 are modeled, though this extrapolation will 
likely breakdown where n e // approaches -3. Simulations that 
model higher particle numbers (leading to higher redshifts) , 
are needed to extend the known range of the mass function of 
dark matter haloes. The accuracy of the S-T function for low 
mass haloes out to high redshifts has important implications 
for a number of astrophysical problems. Evolution of the 
mass and luminosity functions down to dwarf scales permits 
comparison with surveys, providing important cosmological 
tests, and allows calculations of merger histories and star 
formation histories of galaxies, groups, and clusters. Our re- 
sults verify that the S-T function is accurate over the entire 
evolutionary range (for which progenitors or descendants are 
observable) of lyman break galaxies and groups of galaxies. 
Assuming a weak dependence on n e //, the redshift invari- 
ance of the mass function implies that extrapolation of mass 
functions should be reliable for combinations of masses and 
redshifts that cannot presently be simulated, as long as only 
values of lncr -1 that have been verified by simulation are 
considered. The number densities of low mass (< 1O 1O M ) 
haloes at high redshift (z ~ 10), needed for studies of reion- 
ization, or galaxy formation (e.g. Haiman 2002), should be 
well described by the S-T function since although they lie 
below our mass range, they are within our range of simu- 
lated ln<r~ . For such extrapolation to be accurate down to 
indefinitely small masses, all mass would have to be in dark 
matter haloes of some mass or else low mass haloes would 
be overpredicted. 

Extrapolation of mass functions to large values of lncr -1 
that have yet to be verified by simulations, however, are 
likely to be significantly in error, as suggested by the trend 
of increasing overprediction by the S-T function for high 
values of lncr -1 . Though our results only reach ~4cr density 
peaks, there is a trend for the S-T function to increasingly 
overpredict the mass function for increasing a beginning at 
~3<t. The discrepancy with the S-T function for rare objects 
has significant implications for studies that make use of such 
rare objects as a cosmological probe. For example, the num- 
ber density of high redshift (z ~ 6) QSOs, which are thought 
to be hosted by haloes at 5a peaks in the fluctuation field 
(Haiman & Loeb 2001; Fan et al. 2001), are likely to be 
overpredicted by at least a factor of 50%. Some uncertainty 
is also introduced for studies employing the abundance of 
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the highest redshift clusters to probe cosmological param- 
eters (e.g. Robinson, Gawiser, & Silk 2000 and references 
therein). 

5.1 Summary 

In summary, we have utilized high resolution simulations to 
derive the mass function of dark matter haloes over an ex- 
tended range in both mass and redshift over previous work. 
We find that the S-T mass function holds up exceptionally 
over more than 10 orders of magnitude of M/M*. For -1.7 
< In a' 1 < 0.5, the S-T mass function is an excellent fit 
to our data, but begins to overpredict haloes for \n(a~ 1 ) ^ 
0.5, or M/M» ~ 10 B in our volume, reaching a ~50% dis- 
crepancy by ln(cr _1 ) ~ 0.9, corresponding to M/M* ~ 10 9 
at z~15 in our volume. We offer an empirical adjustment 
for the high ln(<r -1 ) portion of the S-T mass function. Our 
results confirm the redshift invariance of the mass function. 
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